home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / GAP.TST / GAPCHISQ.C < prev    next >
C/C++ Source or Header  |  1996-02-12  |  3KB  |  114 lines
  1. /* ============ */
  2. /* gapchisq.c    */
  3. /* ============ */
  4. #include <defcodes.h>
  5. #include <miscdefs.h>
  6. #include <gapdefs.h>
  7.  
  8. static long GapCtrs[MAX_GAP_LEN+1];
  9. /* ==================================================================== */
  10. /* TabulateGapSet - Tabulates Complete Set of Gaps            */
  11. /* ==================================================================== */
  12. /* -------------------------------------------------------------------- */
  13. /* This code reflects the algorithm given by D. E. Knuth, The Art of    */
  14. /* Computer Programming, Vol. 2, Seminumerical Algorithms, Second Ed.,    */
  15. /* Addison-Wesley, Reading (1981), pp. 60-61, Algorithm G.        */ 
  16. /* -------------------------------------------------------------------- */
  17. static void
  18. TabulateGapSet(GAP_DATA_STRU * GapData)
  19. {
  20.     long    GapCt = 0;
  21.     int     GapSize;
  22.  
  23.     /* ----------------------------------- */
  24.     /* Clear Gap Counters & Variate Counts */
  25.     /* ----------------------------------- */
  26.     memset(GapCtrs, 0, sizeof(GapCtrs));
  27.     GapData->TotNumGen    = 0;
  28.     GapData->ActMaxPerGap = 0;
  29.  
  30.     do
  31.     {
  32.     /* ---------------- */
  33.     /* Generate New Gap */
  34.     /* ---------------- */
  35.     GapSize = GenGapData(GapData);
  36.  
  37.     if (GapSize > GapData->ActMaxPerGap)
  38.     {
  39.         GapData->ActMaxPerGap = GapSize;
  40.     }
  41.  
  42.     /* ------------------------------ */
  43.     /* Bump Number Variates Generated */
  44.     /* ------------------------------ */
  45.     GapData->TotNumGen += GapSize + 1;
  46.  
  47.     /* --------------------- */
  48.     /* If No Gap Found, Punt */
  49.     /* --------------------- */
  50.     if (!GapData->CallStatusOK)
  51.     {
  52.         break;
  53.     }
  54.  
  55.     if (GapSize > GapData->MaxGapLen)
  56.     {
  57.         GapSize = GapData->MaxGapLen;
  58.     }
  59.  
  60.     ++GapCtrs[GapSize];
  61.     ++GapCt;
  62.     }
  63.     while (GapCt < GapData->NumGaps);
  64. }
  65. /* ==================================================================== */
  66. /* CalcGapChiSq - Gets Chi-Square Statistic for Gap Test        */
  67. /* ==================================================================== */
  68. void
  69. CalcGapChiSq(GAP_DATA_STRU * GapData)
  70. {
  71.     TabulateGapSet(GapData);
  72.  
  73.     if (GapData->CallStatusOK)
  74.     {
  75.     int    j;
  76.  
  77.     /* ----------------------------------- */
  78.     /* Lump Data in Categories as Required */
  79.     /* ----------------------------------- */
  80.     for (j = 0; j < GapData->NotEnough; ++j)
  81.     {
  82.         GapCtrs[j + 1] += GapCtrs[j];
  83.     }
  84.  
  85.     GapData->GapChiSq = 0;
  86.  
  87.     for (j = 0; j < GapData->NumCategories; ++j)
  88.     {
  89.         if (GapCtrs[j])
  90.         {
  91.         GapData->GapChiSq +=
  92.             SQR((double) GapCtrs[j]) / GapData->CellExpect[j];
  93.         }
  94.         P(printf("\tGapChiSq = %.15e\n", GapData->GapChiSq));
  95.     }
  96.  
  97.     GapData->GapChiSq -= (double)GapData->NumGaps;
  98.  
  99. #    if defined(PRT_GAP_CTRS)
  100.             {
  101.         static x = 3;
  102.         if (x >= 1)
  103.         {
  104.             for (j = 0; j < GapData->NumCategories; ++j)
  105.             {
  106.                 (printf("GapCtrs[%2d] = %ld\n", j, GapCtrs[j]));
  107.             }
  108.             --x;
  109.         }
  110.         }
  111. #    endif
  112.     }
  113. }
  114.